#include "../PhysicalField/PhysicalField2D.h"

#include <cmath>

#include <iostream>
#include <sstream>

#include "../Input/PicParams.h"
#include "../Field/Field2D.h"

#include "../Parallel/ParallelMPI.h"
#include "../Parallel/ParallelMPI2D.h"


using namespace std;

//---------------------------------------------------------------------------------------------------------------------
//Constructor for Electromagn2D
//---------------------------------------------------------------------------------------------------------------------
PhysicalField2D::PhysicalField2D(PicParams* params_, ParallelMPI* pmpi_) :
    PhysicalField(params_, pmpi_)
{
    ParallelMPI2D* pmpi2D = static_cast<ParallelMPI2D*>(pmpi);

    dx       = params->cell_length[0];
    dt_ov_dx = params->dt / dx;
    dx_ov_dt = 1.0 / dt_ov_dx;

    dy       = params->cell_length[1];
    dt_ov_dy = params->dt / dy;
    dy_ov_dt = 1.0 / dt_ov_dy;

    Ex_  = new Field2D<double>(params->dims_field, "Ex" );
    Ey_  = new Field2D<double>(params->dims_field, "Ey" );
    Ez_  = new Field2D<double>(params->dims_field, "Ez" );
    Bx_  = new Field2D<double>(params->dims_field, "Bx" );
    By_  = new Field2D<double>(params->dims_field, "By" );
    Bz_  = new Field2D<double>(params->dims_field, "Bz" );

    rho_  = new Field2D<double>(params->dims_field, "Rho" );
    rho_avg  = new Field2D<double>(params->dims_field, "Rho_avg" );

    //Allocation of time-averaged EM fields
    phi_avg = new Field2D<double>(params->dims_field, "Phi_avg");
    Ex_avg  = new Field2D<double>(params->dims_field, "Ex_avg");
    Ey_avg  = new Field2D<double>(params->dims_field, "Ey_avg");
    Ez_avg  = new Field2D<double>(params->dims_field, "Ez_avg");
    Bx_avg  = new Field2D<double>(params->dims_field, "Bx_avg");
    By_avg  = new Field2D<double>(params->dims_field, "By_avg");
    Bz_avg  = new Field2D<double>(params->dims_field, "Bz_avg");



    rho_global = new Field2D<double>(params->dims_field_global, "Rho_global");
    phi_global = new Field2D<double>(params->dims_field_global, "Phi_global");
    Ex_global = new Field2D<double>(params->dims_field_global, "Ex_global");
    Ey_global = new Field2D<double>(params->dims_field_global, "Ey_global");
    Ez_global = new Field2D<double>(params->dims_field_global, "Ez_global");


    rho_global_avg = new Field2D<double>(params->dims_field_global, "Rho_global_avg");
    phi_global_avg = new Field2D<double>(params->dims_field_global, "Phi_global_avg");
    Ex_global_avg  = new Field2D<double>(params->dims_field_global, "Ex_global_avg");
    Ey_global_avg  = new Field2D<double>(params->dims_field_global, "Ey_global_avg");

    rho_global->set_value(0.0);
    phi_global->set_value(0.0);

    Ex_->set_value(0.0);
    Ey_->set_value(0.0);
    Ez_->set_value(0.0);
    Bx_->set_value(params->externB[0]);
    By_->set_value(params->externB[1]);
    Bz_->set_value(params->externB[2]);
    rho_->set_value(0.0);



    //Allocation of the time-averaged EM fields
    Ex_avg  = new Field2D<double>(params->dims_field, "Ex_avg" );
    Ey_avg  = new Field2D<double>(params->dims_field, "Ey_avg" );
    Ez_avg  = new Field2D<double>(params->dims_field, "Ez_avg" );
    Bx_avg  = new Field2D<double>(params->dims_field, "Bx_avg" );
    By_avg  = new Field2D<double>(params->dims_field, "By_avg" );
    Bz_avg  = new Field2D<double>(params->dims_field, "Bz_avg" );

    //Charge currents currents and density for each species
    for (int ispec=0; ispec<n_species; ispec++)
    {
        rho_s[ispec]            = new Field2D<double>(params->dims_field, ("Rho_"+params->species_param[ispec].species_name).c_str());
        rho_s_avg[ispec]        = new Field2D<double>(params->dims_field, ("Rho_"+params->species_param[ispec].species_name+"_avg").c_str());

        rho_s_global[ispec]     = new Field2D<double>(params->dims_field_global, ("Rho_global_"+params->species_param[ispec].species_name).c_str());
        rho_s_global_avg[ispec] = new Field2D<double>(params->dims_field_global, ("Rho_global_"+params->species_param[ispec].species_name+"_avg").c_str());

        Vx_s[ispec]             = new Field2D<double>(params->dims_field, ("Vx_"+params->species_param[ispec].species_name).c_str());
        Vx_s_avg[ispec]         = new Field2D<double>(params->dims_field, ("Vx_"+params->species_param[ispec].species_name+"_avg").c_str());
        Vx_s_global[ispec]      = new Field2D<double>(params->dims_field_global, ("Vx_global_"+params->species_param[ispec].species_name).c_str());
        Vx_s_global_avg[ispec]  = new Field2D<double>(params->dims_field_global, ("Vx_global_"+params->species_param[ispec].species_name+"_avg").c_str());

        Vy_s[ispec]             = new Field2D<double>(params->dims_field, ("Vy_"+params->species_param[ispec].species_name).c_str());
        Vy_s_avg[ispec]         = new Field2D<double>(params->dims_field, ("Vy_"+params->species_param[ispec].species_name+"_avg").c_str());
        Vy_s_global[ispec]      = new Field2D<double>(params->dims_field_global, ("Vy_global_"+params->species_param[ispec].species_name).c_str());
        Vy_s_global_avg[ispec]  = new Field2D<double>(params->dims_field_global, ("Vy_global_"+params->species_param[ispec].species_name+"_avg").c_str());

        Vz_s[ispec]             = new Field2D<double>(params->dims_field, ("Vz_"+params->species_param[ispec].species_name).c_str());
        Vz_s_avg[ispec]         = new Field2D<double>(params->dims_field, ("Vz_"+params->species_param[ispec].species_name+"_avg").c_str());
        Vz_s_global[ispec]      = new Field2D<double>(params->dims_field_global, ("Vz_global_"+params->species_param[ispec].species_name).c_str());
        Vz_s_global_avg[ispec]  = new Field2D<double>(params->dims_field_global, ("Vz_global_"+params->species_param[ispec].species_name+"_avg").c_str());


        Vp_s[ispec]             = new Field2D<double>(params->dims_field, ("Vparallel_"+params->species_param[ispec].species_name).c_str());
        Vp_s_avg[ispec]         = new Field2D<double>(params->dims_field, ("Vparallel_"+params->species_param[ispec].species_name+"_avg").c_str());
        Vp_s_global[ispec]      = new Field2D<double>(params->dims_field_global, ("Vparallel_global_"+params->species_param[ispec].species_name).c_str());
        Vp_s_global_avg[ispec]  = new Field2D<double>(params->dims_field_global, ("Vparallel_global_"+params->species_param[ispec].species_name+"_avg").c_str());

        T_s[ispec]              = new Field2D<double>(params->dims_field, ("T_"+params->species_param[ispec].species_name).c_str());
        T_s_avg[ispec]          = new Field2D<double>(params->dims_field, ("T_"+params->species_param[ispec].species_name+"_avg").c_str());
        T_s_global[ispec]       = new Field2D<double>(params->dims_field_global, ("T_global_"+params->species_param[ispec].species_name).c_str());
        T_s_global_avg[ispec]   = new Field2D<double>(params->dims_field_global, ("T_global_"+params->species_param[ispec].species_name+"_avg").c_str());
    }

}



PhysicalField2D::~PhysicalField2D()
{

}



void PhysicalField2D::incrementAvgFields(int itime)
{
    //log_info<<"incrementAvgFields 2D";
    if( (itime - 1) % params->timestep_io == 0 )
    {
        rho_global_avg->set_value(0.0);
        phi_global_avg->set_value(0.0);
        Ex_global_avg->set_value(0.0);
        Ey_global_avg->set_value(0.0);

        for(int ispec=0; ispec<n_species; ispec++)
        {
            rho_s_avg[ispec]->set_value(0.0);
        }
    }

    //Calculate the sum values for global rho phi Ex and Ey
    if( (itime % params->timestep_io) > (params->timestep_io - params->timestep_average) || (itime % params->timestep_io) == 0 )
    {
        for(int i = 0 ; i < dims_totalglobal ; i++)
        {
            (*rho_global_avg)(i) += (*rho_global)(i);
            (*phi_global_avg)(i) += (*phi_global)(i);
            (*Ex_global_avg)(i)  += (*Ex_global)(i);
            (*Ey_global_avg)(i)  += (*Ey_global)(i);
        }

        //Calculate the sum values for density of each species
        for(int ispec = 0; ispec < n_species; ispec++)
        {
            //all fields are defined on the primal grid
            for(int i = 0 ; i <dims_total; i++)
            {
                (*rho_s_avg[ispec])(i) += (*rho_s[ispec])(i);
            }
        }
    }

    //calculate the averaged values
    if( itime % params->timestep_io == 0 )
    {
        for(int i=0 ; i<dims_totalglobal ; i++)
        {
            (*rho_global_avg)(i) /= params->timestep_average;
            (*phi_global_avg)(i) /= params->timestep_average;
            (*Ex_global_avg)(i)  /= params->timestep_average;
            (*Ey_global_avg)(i)  /= params->timestep_average;
        }
        for(int ispec = 0; ispec < n_species; ispec++)
        {
            for(int i = 0 ; i < dims_total ; i++)
            {
                (*rho_s_avg[ispec])(i) /= params->timestep_average;
            }
        }
    }

}
